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Abstract Simulations of monodisperse and polydis- 
perse (112(A) = 0.13 ±0.002) 2D foam samples undergo- 
ing simple shear are performed using the 2D Viscous 
Froth (VF) Model. These simulations clearly demon- 
strate shear localisation. The dependence of localisation 
length on the product XV (shearing velocity V times 
external wall friction coefficient A) is examined and is 
shown to agree qualitatively with other published ex- 
perimental data. A wide range of localisation lengths is 
found at low XV, an effect which is attributed to the 
existence of distinct yield and limit stresses. The gen- 
eral Continuum Model is extended to incorporate such 
an effect and its parameters are subsequently related to 
those of the VF Model. A Herschel-Bulkley exponent of 
a = 0.3 is shown to accurately describe the observed 
behaviour. The localisation length is found to be inde- 
pendent of XV for monodisperse foam samples. 
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1 Introduction 

Foam is defined as a two-phase system in which a dis- 
persed phase of gas is enclosed by a continuous phase 
of liquid [25]. In aqueous foams, the dispersed phase is 
typically air, and the liquid phase water with an added 
surfactant. To simplify the task of understanding the 
rheology of these systems, it has become popular to con- 
centrate on two-dimensional (2D) foams, which consist 
of a single planar layer of bubbles. It is the rheology 
of these systems which is under scrutiny in this paper, 
using the 2D VF Model. 

We endeavour to understand the mechanisms which 
cause localisation of shear at a moving boundary, as re- 
ported in a number of recent experiments (see below). 
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This behaviour is observed in the VF simulations that 
we will discuss. In the literature, the extent of this shear 
localisation effect is measured by a localisation length. In 
this paper, two different definitions of localisation length 
will be employed. 

Simulation results for polydisperse froths subjected 
to simple shear will be presented, where the localisa- 
tion length is found to vary with XV (viscous drag times 
driving velocity; why this is the important parameter 
to consider is explained in Sec. [4]). For low XV, a wide 
range of localisation lengths is found; see Sec. |3J This 
effect is attributed to the existence of distinct yield and 
limit stresses. We give qualitative evidence of this in Sec. 
[5j where we extend the general Continuum Model to in- 
corporate such an effect. In Sec. [6] we proceed to relate 
parameters of the VF Model to the Continuum Model 
which is shown to accurately predict the observed be- 
haviour for a Herschcl-Bulkley exponent of a = 0.3. The 
localisation length is found to be independent of XV in 
the monodisperse case. 

In experiments with 2D foams, often a single layer of 
bubbles is confined between two narrowly spaced glass 
plates (a Hele-Shaw cell). There are also other types of 
quasi-2D systems, such as the Bragg raft [3], where a sin- 
gle layer of bubbles floats on a liquid pool, and the con- 
fined bubble raft, where a Bragg raft is trapped under- 
neath a glass plate. The most important distinction be- 
tween these experimental realisations is concerned with 
the presence of viscous drag. When a foam is in contact 
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with one or two confining plates, there is a drag force 
associated with any movement of the foam relative to 
the plate(s). As we shall see, the drag force in the VF 
Model (see eq. [2| plays a role in the localisation of flow 
in our foam samples. 




Fig. 1 A Tl neighbour-swapping event triggered by apply- 
ing shear, a) Initial configuration, b) A and B lose their com- 
mon edge, creating an unstable four-fold vertex point, c) a 
new edge is created between C and D and d) final configura- 
tion. Data taken from a VF simulation. 

When a 2D foam is subjected to an applied shear 
stress, after an initial transient, it yields and begins to 
flow. The foam yields locally when the yield stress is 
reached. A more detailed description of the stress-strain 
relation will be required when interpreting the simula- 
tion results presented in this paper; see Sec. [5] At the 
local level, yielding is due to plastic events, i.e. Tl topo- 
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logical changes of the foam structure (see Fig. [IJ . Two 
neighbouring bubbles (A and B) lose a common edge 
which is subsequently gained by two proximate bubbles 
(C and D), which become neighbours. We describe the 
incorporation of these topological changes into the VF 
Model in Sec. [2] 

When flow is concentrated in one region and not in 
another, the flow is said to have localised. Debregeas et 
al. [8] were the first to report definitive evidence of shear 
localisation in 2D aqueous foams. Their experiments ex- 
hibit shear localisation next to a moving boundary in a 
Couette geometry, with an exponential decay in the mea- 
sured foam velocity profiles. Similar results have been re- 
ported by Wang et al. [25] and Krishan and Dennin [23] 
for straight and circular geometries, respectively. These 
results have been interpreted within the framework of 
the Continuum Model [6l[T4l[T5] where shear localisation 
is attributed to the presence of a drag force. This notion 
is further supported by the work of |12j which studies the 
effects of drag forces at high shear rates using the VF 
Model. However there are also quasi-static simulations 
showing localisation (discussed below) in which there is 
no such wall drag. This suggests that there is more than 
one mechanism that may lead to shear localisation. We 
will return to this point in Sec. [5] 

Experimental work by Katgert et al. I 9, [3D] on the 
shearing of bidisperse foams in a Hele-Shaw cell (straight 
geometry, that is, simple shear) shows Herschcl-Bulkley 
behaviour (discussed below) and supplies further evi- 
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dence of shear localisation in 2D foams. In this case, 
however, the velocity profiles are not exponential. Such 
non-exponential velocity profiles, together with their ve- 
locity dependence can be obtained from an extension or 
generalisation of the original Continuum Model [35] [30] ■ 
Furthermore, these experiments show that the locali- 
sation length decreases as the velocity of the moving 
boundary increases. In this paper, we show velocity pro- 
files from our VF simulations which exhibit qualitatively 
similar behaviour. 

Of particular current interest in these types of exper- 
iments is the dependence of (local) shear stress on shear 
rate. This effect is captured by the Herschel-Bulkley con- 
stitutive relation, 

cr = (j y + c v e a (1) 

where a is stress, a y is the yield stress, the coefficient 
c v is the so-called consistency, e is strain rate and a is 
the HB exponent. Katgert and co-workers report a = 
0.36. They also note that in the monodisperse case (i.e. 
bubbles of equal size) , the localisation length is found to 
be independent of shear rate. We too find this to be the 
case in our simulations; see Sec. [4] 

Shear localisation has been studied computationally 
using other microscopic (bubble scale) models. Quasi- 
static models, as explored by [5] [T3J HU GH] might shed 
light on behaviour at very low strain rates. Results re- 
ported by Kabla [121 HE] show localisation next to the 
boundaries in quasi-static shearing simulations (^{A) ~ 
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0.06 using the definition eq. 131), consistent with experi- 
mental observation. In these simulations of simple shear, 
where there is no wall drag, either wall may be regarded 
as the one that moves. Recent results by Wyn [33] sug- 
gest that as the second moment of the bubble area dis- 
tribution 112(A) is increased in such simulations to val- 
ues approaching 0.2 or higher, shear-banding can occur 
in regions away from the moving boundary. The width 
of these shear bands has a square-root dependence on 
H2 (A) . This type of behaviour has yet to be observed in 
experiments. In this paper, we only examine the rheol- 
ogy of foam samples of disorder 112(A) = 0.13 but will 
probe higher values of disorder in future work in order 
to search for similar effects. 

With the aid of the Surface Evolver software [J], 
quasi-static simulations are increasingly easy to imple- 
ment and, with modern computing, are certainly fast. 
But are they suitable for rheology? In quasi-statics, the 
foam is relaxed to equilibrium at each step. There is 
therefore no relevant time scale present and so no con- 
cept of shear-rate. It makes no sense to consider Herschel- 
Bulkley type relations or to discuss the dependence of 
localisation on boundary velocity. 

What then, are the alternatives to quasi-static sim- 
ulations? Bubble models [5], where a foam is modeled 
as a collection of interacting disks appear to represent 
at some level the dynamics of 2D foams. Langlois et 
al. [H] report a Hcrschcl-Bulkley exponent of a = 0.54 
(112(A) ~ 0.03). In addition, shear localisation is ob- 
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served when wall drag is present. For dry foams though, 
where low liquid fraction causes bubbles to become more 
polygonal in shape, this approach is no longer accurate 

M- 

In this paper, we adopt the 2D Viscous Froth (VF) 
Model [H] [21] as a more realistic model for dry 2D 
foam dynamics. We have performed an extensive study 
of shear localisation with the VF Model in a straight ge- 
ometry which shows realistic dynamics and a rich variety 
of behaviour, particularly at low AV. 

For a summary of the experimental and theoretical 
work presented in this section, see |27) . 

2 The 2D Viscous Froth Model and its 
implementation 

The model describes the motion of a soap film in the 
2D systems described above, with wall drag [21]. In the 
present case, bubble areas are kept constant. The foam is 
considered to be sufficiently dry (liquid fraction less than 
0.01) so that a soap film may be accurately described by 
a curved line and the junctions are represented by points. 
In the present simulations, a soap film is approximated 
as a system of connected straight line segments. The 
motion of a point s joining these segments is given via 
the equation 

Xv ± (s) = AP-jK(s) (2) 
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Fig. 2 A diagram illustrating the various forces involved in 
film motion with the two-dimensional VF Model. B and B' 
indicate the two bubbles the central soap film is separating. 
Note that this film is in contact with a surface (in the plane 
of the page) which results in a drag force when the film is 
moving. 

where A is the wall drag coefficient, v ± (s) is the velocity 
of a point s in the direction of the normal vector N(s) 
to the soap film, AP represents relative pressure differ- 
ences between neighbouring cells, 7 is a constant surface 
tension force (in 2D), and K(s) is local film curvature 
calculated from the relative positions of adjacent (dis- 
crete) film segment points. See Fig. [2] for an illustration 
of the forces involved. Setting A = in eq. [2] recovers the 
Young-Laplace law, corresponding to soap films that are 
arcs of circles. 

Throughout the implementation of the model, film 
segments adjacent to the three-fold vertex points are 
held an angle of ~ radians relative to each other, in 
accordance with Plateau's rules for a soap froth. Dc- 
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tails on the numerics of this calculation are best found 
in [TT]. It should be noted that for high rates of strain, 
one would expect surface tensions in the soap films to 
vary to the point that this equilibrium condition would 
no longer apply (for example, because of the M arangoni 
effect). At least for lower rates of strain, the ^ rule is 
reasonable. 

The VF model may be conveniently incorporated into 
a Surface Evolver 0] script (as pioneered by Cox [7]), 
thus allowing for the use of various SE features. The 
procedure for performing (Tl) topological changes in the 
Surface Evolver is as follows, and is illustrated in Fig. [I] 
The distance along the film between neighbouring three- 
fold vertex points is calculated at each timestep. When 
this film length becomes smaller than a predefined criti- 
cal cut-off length, l c then it is deleted using the Evolver's 
'edgeweed' command. A four-fold vertex is temporarily 
formed to maintain the topology of neighbouring cells 
(see Fig. [ljb)). The Evolver's 'pop' command is then 
employed, which scans the foam for vertices which do 
not have a legal topology and replaces the four-fold ver- 
tex with a proper local topology. This results in a new 
film of effectively negligible length oriented in the per- 
pendicular direction to the old film (see Fig.Qc)). 

Further details on the implementation of the VF Model 
can be found in the papers by Kern et al. [21] and Green 
et al. [TT] . 
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3 Sample Creation 

A semi-periodic monodisperse sample is created using 
the standard method outlined in the Surface Evolver 
documentation (which is supplied with the software pack- 
age). Disordered semi-periodic samples are created by 
the following process (illustrated in Fig. |3|. 

Points are placed at random in the unit cell using 
a uniform distribution to determine the x and y posi- 
tions; see Fig. [3]ja). New points are added to the box if 
they are more than a predefined minimum distance r m i n 
from any other point. The process is continued until the 
desired number of points have been successfully placed. 
A lower value of r m i n results in more polydisperse sam- 
ples. These points are translated to boxes to the left 
and right, and reflected (see Fig. [3]ja)) through the lines 
y = and y = 1 to boxes above and below, based on 
the method of De Fabritiis and Coveney [TU] (as indi- 
cated by the background trian gles in Fig. [jjja)). With 
the software package Qhull pQ, the Voronoi Diagram (a 
particular way of tessellating the plane into regions of 
convex polygons) of these points is calculated; see Fig. 
[3]jb). This is then passed into the Surface Evolver. The 
box in the centre is isolated (see Fig. [3^c)) and keeping 
the areas of each of the cells fixed, the Surface Evolver 
performs line minimization on the structure. The result- 
ing structure is our final two-dimensional half-periodic 
(i.e. periodic in the x-direction only) foam data file; see 
Fig. ^d) . Of interest here (as a result of the reflection) 
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is that the straight line boundaries at y = and y = 1 
naturally occur as a result of this process. 

4 Simulation Details and Results 

Using the above methods for foam sample creation, wc 
create one monodisperse foam sample and five foam sam- 
ples of polydispersity p, 2 (A) = 0.13 ± 0.002, where the 
measure of polydispersity is defined as the second mo- 
ment of the area distribution, 

V2(A)=(l-jJ. (3) 

Here A denotes the area of a bubble, and A the mean 
bubble area. 

The foam samples consist of Nf, = 100 bubbles in 
a square unit cell of area 1; it is too computationally 
expensive to run larger samples in a VF simulation. In 
our dimensionless simulation units, our system size L = 
1 and mean bubble area A = 0.01. We define a new 
length scale A 1 / 2 , the square root of the mean bubble 
area. In these new units, L = 10 A 1 / 2 . The width Wj of 
one layer of bubbles in our square sample is given by 

Wi = L/^N~ b = A 1 ' 2 (4) 

We proceed to move the top boundary in the posi- 
tive x-direction with velocity V by incrementally moving 
vertices at y = 1 a distance Vdt per timestep dt. The 
VF algorithm, as outlined in Sec. [2] is used to determine 
the dynamics of the foam during each timestep. Typical 
values for the displacement of the shearing boundary per 




Fig. 3 The creation of a semi-periodic polydisperse two-dimensional foam sample as required for our simple shear simulations, 
a) Points are placed in the central unit cell and translated/reflected into adjacent boxes (as indicated by the background 
triangles), b) The Voronoi Diagram of these points is calculated, c) The central area is isolated and made into a half-periodic 
(in x-direction) data-file, d) Keeping cell areas constant, the Surface Evolver performs line minimisation on the structure. We 
refer to this shown structure as Sample 1 later in the text. 

timestep are in the range lCT 6 ^ 1 / 2 < Vdt < lCT 3 ,? 1 / 2 

(depending on what values of V and A are used) . No-slip 

boundary conditions are maintained by fixing vertices ly- 

Multiple simulations are run for different values of 

ing on the boundaries while the VF algorithm is being 

XV (wall drag coefficient times boundary velocity) with 

implemented. Our boundary conditions are thus 

a fixed value of surface tension 7. To see why this is 
the appropriate parameter to look at, consider again the 
equation of motion for the VF Model, as given by eq. [2] 




(5) 
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By setting v — Vv , where V is the boundary velocity 
and v 1 - is our rescaled dimensionless velocity, we can 
rewrite our equation of motion as 

(AF)u- L (s) = AP - jK(s) (6) 

It is clear that, given any initial state configuration, its 
development in time is determined by XV. Furthermore, 
we see evidence of this XV dependence if we rewrite the 
Herschel-Bulkley relation (see eq. [I]) in terms of our VF 
parameters. As stress in 2D has dimensions of force per 
length, on dimensional grounds, we see that 

a = <7 y + c v 7 1 ~ a A a ~ 1 / 2 L~ a (XV) a (7) 

where the 2D surface tension 7 has dimensions of force, 
XV has dimensions of force per length and c v is a dimen- 
sionless parameter of order unity which may be related 
to (12(A). In this derivation, we define the strain rate 
term of eq. [T] as the nominal shear rate of the system, 
e = V/L. 

To calculate flow profiles, bubble centre positions are 
determined. We subsequently divide our foam into bins 
of width Wi and calculate the average velocity of bubbles 
centres in each bin over time. A sketch of our simulation 
setup is illustrated in Fig. [4] (polydisperse sample) . 

Fig. [5] shows examples of averaged steady state ve- 
locity profiles. We say that a simulation has reached 
a steady state once there is no longer any appreciable 
change in our velocity profile in time. Typically, we av- 
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Fig. 4 (a) A polydisperse foam consisting of 100 bubbles 
(Sample 5) in equilibrium. L denotes our system size, (b) 
The same foam being sheared at a velocity V. 

erage our steady state velocity profiles over the range 
1 < e < 10, where the imposed strain e is defined as 
e = Ax I L and Ax is equal to the total displacement of 
the moving boundary. Note that there is a clear change in 
the flow profiles as we vary XV. We find that localisation 
occurs close to the moving boundary in all but two of our 
simulations. (In one of these cases, for XV = 0.01 ^A 1 / 2 , 
localisation switches to the stationary boundary, while in 
the second case, for XV — 0.005 7A 1 / 2 , a shear band oc- 
curs in the centre of the sample away from either bound- 
ary (data not shown).) Localisation of flow can also be 
made visible by plotting the positions at which Tl topo- 
logical changes occur in our samples, as done by [33]. An 
example is shown in Fig. [6j 

At this stage, it is unclear what the form of the ve- 
locity profiles is. We have attempted to use exponential 
fits and fits from the general Continuum Model [28] in 
the data fitting process but this approach does not yield 
consistently good fits to our velocity profiles which are 
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= 0.01 y/A 1 ' 2 


++XV 


= 0.025 y/A" 2 




= 0.05 y/A 1 ' 2 


mm xv 


= 0.1 y/A" 2 




clearly quite noisy, presumably due to the small system 
size. To obtain a measure of the width of the flowing 
region from these noisy profiles, we use the following 
definition of localisation length, denoted by U n t [2Z] 



1 

V 



v(y)dy . 



(8) 



2 4 6 

— 1/2 

Distance from moving boundary [A ] 



Fig. 5 Examples of velocity profiles for different values of 
XV. Profiles shown are for Sample 1. The corresponding lo- 
calisation lengths for these profiles are denoted by filled cir- 
cles in Fig. [7] 
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Fig. 6 The location of Tl topological changes as a function 
of strain for a polydisperse sample. After an initial transient, 
where Tls happen everywhere in the foam, the flow localises 
and Tls are found to occur mostly next to the moving bound- 
ary (at y = 10 A 1 / 2 ). Shown data is for Sample 1 where 
XV = 0.05 7/A 1/2 



Fig. 7 Localisation lengths for a range of XV. Each symbol 
represents one simulation run. For low XV, there is a wide 
range of lengths, while for high XV, only the first layer of 
bubbles flows. Filled symbols indicate simulations where V 
is fixed and A is varied. Open symbols indicate simulations 
where A is fixed and V is varied. 

This integral, which has the required dimensions of 
length, is calculated numerically for each of our velocity 
profiles using the Trapezoidal Rule. Fig.[7]shows a varia- 
tion of localisation length with XV. Note that for low XV 
we find large scatter in the localisation lengths, however, 
this scatter decreases as XV is increased. For high XV, 
the length converges towards the minimum localisation 
length, l min = A 1 / 2 , the width of one bubble layer (see 
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eq. [1| . This is because the first layer of bubbles always 
flows. 

The ratio of the intrinsic timescale of the VF Model 
to the external timescale (as imposed by the nominal 
shear rate e = V/L), otherwise known as the Deborah 
number _D e , is given by 



as defined in [3T]. A small Deborah number (D e <C 1) 
indicates that the foam has enough time available to re- 
equilibrate, even as the applied shear attempts to bring 
the foam out of equilibrium (and vice versa for large D e ). 
In our simulations, 0.001 < D e < 0.03. We therefore 
conclude that we are close to the quasi-static regime in 
all of the discussed simulations. 

While not shown here, similar simulation runs have 
been performed for a monodispcrsc foam (^(A) = 0). 
The localisation length is found to be independent of XV 
and is determined to be I — A 1 / 2 (the same as l m in in 
our polydisperse simulations). 

Our simulation results are broadly consistent with 
the findings of Katgert et al. [THl HO] , where the localisa- 
tion length is found to decrease with increasing V, and 
rate independence of localisation length is found in the 
monodisperse case. 

We wish to gain an understanding of these VF sim- 
ulation results by attempting to capture the observed 
behaviour by a continuum model. Such a model must 
include a constitutive relation which relates the local 
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(wall) drag force of the VF model to an averaged drag 
force in the continuum description. However, this would 
not be enough to explain the observed simulation results, 
as (according to the general Continuum Model [28]) it 
would result in a zero localisation length in all cases. 
Therefore, results suggest that internal dissipation (rep- 
resented by the shear rate term in the HB relation; see 
cq. [T]) should also be included, although it is not clear 
how this dissipation arises in the VF simulations. We 
will also appeal to the idea of the existence of a stress 
overshoot in order to explain the variation of localisation 
length at low XV. 

5 The Continuum Model 

Up until now, we have discussed microscopic (bubble 
scale) models of 2D foam rheology. An alternative way 
of describing a foam is to treat it as a continuum. The 
generalised Continuum Model [35] [50] (for steady shear) 
combines the Herschel-Bulklcy constitutive relation (see 
eq. [l]) with the following expression for the variation of 
(wall) drag force Fd per unit area as a function of local 
velocity v 

F d = -c d v b , (10) 

where Cd is the drag coefficient and b is the drag exponent 
(the Bretherton law gives b = | [5]). These two expres- 
sions can be related by a force balance, which leads to 
the following differential equation [TH [5U] 
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d 
dy 



dv(y) 



dy 



v(y) b 



(11) 



which can be solved using the boundary conditions v(0) — 
V and v(L) = 0. These are equivalent to the boundary 
conditions imposed in our VF simulations; see eq. [5] (al- 
beit that here the distance y is measured downward from 
the shearing boundary). 

Upon inspection, it is clear that that a velocity profile 
of the following form satisfies eq. and exhibits flow 
localisation: 



v(y) = V(l - y/yoY 



(12) 



where yo and n may be obtained by inserting eq. 12 into 
and equating prefactors and exponents [27] . This 



eq 
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gives 



and 



1 + a fa(l + b)c v V 
a-b \ (1 + a)c d 



1 + a 
a — b 



a—b\ i+a 



(13) 



(14) 



Eq. [12] clearly satisfies the first boundary condition, 
v(0) = V of eq. 11 The second boundary condition is 
satisfied if we take our sample size L — > 00 [27] . This 
approach is valid so long as the size of the sample is 
much greater than the localisation length (L ^> I) which 
may be defined as 



/ = 



V 



dv{0) 



(15) 



an alternative definition to that of the previous section; 
see eq. [8] Inserting eq. [12] into eq. [15] leads to the follow- 
ing expression for localisation length as a function of the 
shearing velocity V, 

\(l + a)c d J 

Its relation to the definition of localisation length l int 
(see eq. |8| is 



km/I = (1 + a)/(l + 2a - b) 



(17) 



as given by [27J. For a < b (as is the case for our VF 
simulations; see Sec. [6]), localisation length therefore de- 
creases as the shearing velocity V is increased. 

The possibility of having a range of localisation lengths 
at low V (as in Fig. [7]) can be accounted for by extending 
the Continuum Model to incorporate what we will refer 
to as a stress overshoot. This we will now proceed to do. 

In a recent paper [28], Weaire et al. introduced the 
idea of distinct yield o~ y and limit cr; stresses as a possi- 
ble mechanism for localisation in the absence of viscous 
drag. An illustration of the typical stress vs strain pic- 
ture is shown in Fig. |8] The constitutive stress relation 
thus becomes 



a = ui + c v e 



(18) 



where 07 denotes the limit stress. When the magnitude of 
the stress overshoot, A = a y —07 is set to zero we recover 
the original Herschel-Bulkley relation (see eq. [T]) . 
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The differential equation given by eq. [TT] may be 
solved numerically, yielding velocity profile solutions of 
the kind we envisage, which are valid if they satisfy the 
inequality given by eq. [f9] As the local strain rate is 
defined as 



Strain 



Fig. 8 An illustration of a stress vs strain relation incorpo- 
rating the idea of distinct yield and limit stresses, denoted 
by G v and ai respectively. The filled circles indicate that the 
foam can co-exist at the same stress at the boundary between 
flowing and non-flowing regions. 

If shear localisation is present in a foam, there exists 
at least one point ys which lies on the boundary between 
flowing and stationary regions. In our VF simulations, 
this corresponds to the point at which the velocity profile 
intercepts the v = axis (see, for example Fig. [5]). At 
this point, the system can co-exist at the same value of 
stress in both static and flowing regions, as indicated by 
the filled dots in Fig. [8} The stress at ub can take on any 
value between ui and o~ y as the foam is sheared. From 



Kv) 



dv(y) 



dy 



(20) 



eq. 18 we see that this leads to the inequality 



< c v e{y B ) a < A 



(19) 



As V is increased, on average we expect the viscous 
stress c v e(y) a between and ys to cause the stress in 
the flowing region to lie closer to a y so that the stress 
overshoot is less evident. However, at low V, the effect 



is obvious (see Fig. 10 1 and may have important effects. 



in this case [27], the quantity e(ys) can be directly mea- 
sured from the calculated velocity profiles, provided they 
intersect the v = axis at some finite value. 

The results of these calculations can be seen in Fig. 
[9] where we have solved the model numerically for the 
values a = 0.5, b — = c v = A = 1. The upper bound 
l + (V) corresponds to where the shear stress ct(?/b) = 07, 
where the viscous stress c v e(yg) a = and the ana- 
lytic solution for localisation length given by eq. [16] The 
lower bound l~ (V) corresponds to where the shear stress 
\vb) — o~y and where the viscous stress c v i{yB) a = A. 
Thus, for a given V, l~(V) < 1{V) < l+(V) gives the 
range of allowed solutions, indicated by the shaded re- 
gion in Fig. [9] 

We note that Fig. [9] is qualitatively similar to Fig. 
[jj with a large range of possible localisation lengths at 
low V and convergent behaviour at high V. Remarkably, 
the model predicts that for low V, the localisation length 
can take any value < I < 00. This prediction, of course, 
holds only in the presence of viscous drag. 

An important question to be answered is how does 
one define the critical velocity V c below which the foam 
can take on a wide range of localisation lengths? If one 
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_ a=0.5 , b=l , c =1, c =1 , A=l 

7 ' d ' V ' 



. g+j fa(l + b)\ 1+6 / 1 

Vr = A .(! + ») f - J — — 

1 + a / \c d {c v )-. 



(22) 



2 3 
V 



Fig. 9 Results from a numerical solution of the Continuum 
Model incorporating the inequality given by eq. |19| For low 
V, there is a large range of possibilities for localisation length, 
while for high V, the range of allowed lengths converges to a 
narrow band. V c is our critical velocity below which a large 



range of localisation lengths is possible; see eq. 22 l + (V) 
and l~ (V) denote the upper and lower bounds to the range 
of allowable solutions, respectively. 

assumes that as V — > 0, the velocity profile becomes 
approximately linear, then e 
isation length. If we are on the lower bound, from eq. 
19 we see that A = c v fy) , or expressing it in a more 
convenient form, 



y, where I is the local- 



l = 



(I) 



V . 



(21) 



We are interested in the point where this line inter- 
sects the upper bound l + (V), which is given by eq. 16 



We solve this pair of simultaneous equations (eq. 16 21 1 



in terms of V and choose to define the point of intersec- 
tion as our critical velocity V c (see Fig. [9]). This yields 



To make a more quantitative comparison between 
continuum theory and the VF results presented in Sec. 
|4j a more detailed study of the relationship between the 
parameters of both models is required. We present such 
a study in the next section. 

In one of the earliest publications on this subject, 
Kabla & Debregeas [IT] attribute shear localisation in 
quasi-statics to what they call 'self amplification'. This 
idea is qualitatively the same as the ideas presented in 
this section. This approach may have the capacity to 
explain other results in the literature, particularly |26j 
where shearing experiments are performed for an ordi- 
nary Bragg raft (where there are no confining plates) 
in a straight geometry. In these experiments (as in our 
VF simulations) variations in the averaged velocity pro- 
files are observed between experiments but averages over 
several experiments converge much better. 

6 Relating the Continuum Model to the Viscous 
Froth Model 

We now proceed to relate the parameters of the (mi- 
croscopic) VF Model and the (macroscopic) Continuum 
Model. This is done using a combination of numerical 
and analytic approximations. 

To demonstrate preliminary evidence of existence of 
the stress overshoot in simulation, we have performed 
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quasi-static calculations using the Surface Evolver (of 
the type mentioned in in Sec. [lj for 29 foam samples of 
disorder fi 2 (A) = 0.13 ± 0.03 with N b = 100 bubbles. 
This effectively sets the viscous stress c v e(y) a to zero, 
thereby allowing us to obtain an accurate estimate of the 
magnitude of the stress overshoot, A. The foam samples 
are created using the process outlined in Sec. [3] The 
shear stress a xy (defined in (55]) is recorded for each 
simulation and subsequently averaged; see Fig. [lO} 

As there is localisation in these simulations (at either 
the moving or stationary boundary) which affects our 
stress measurements, the limit stress o\ reported here 
must be treated as an approximate measurement. The 
value of the yield stress a yi however, is exact, as up to 
a strain of unity, the foam is in the elastic regime and 
the bubble motion has not yet localised. We measure the 
magnitude of the stress overshoot to be A — 0.1 7/j? 1 / 2 , 
which corresponds to a 17% overshoot. In the calculation 
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shown in Fig. 12 a 20% overshoot is used. 



Fig. 10 Averaged shear stress data for 29 foam samples of 
disorder /X2(j4) = 0.13 ± 0.03. The stress overshoot is clearly 
evident. The yield stress a y is taken to be the maximum stress 
value, which occurs at a strain of unity. The limit stress <t; 
is the stress average from a strain of 2 to 10. Calculation 
performed using quasi-static simulations. 

force is defined and the orientation of soap films in the 
foam. 

The drag force per unit area must be proportional to 
the total length of the soap films in that area. For the 
honeycomb, this yields 



The drag force per unit area for the Continuum Model 
is given by eq. [TUJand acts in the direction of shear. We 
wish to relate this to the the drag force of the VF Model, 
Xv 1 , which is a force per length and acts in the normal 
direction to a soap film (see Fig. [2]) . Trivially, the drag 
exponent, 6=1. The numerical prefactor may be 
calculated analytically for a 2D hexagonal honeycomb 
structure, which serves as a reasonable approximation. 
We also take into account the direction in which the drag 



c d oc 



'2V3 



(23) 



In the VF simulations, it is observed that bubbles 
move on average only in the direction of shear. The mag- 
nitude of this 'apparent' velocity is denoted by v app in 



Fig. 11 However, the drag force for the VF model by 
definition points in the direction of the normal to a soap 
film, and so we project v app in this direction (see Fig. 
111)). This results in |v I — I v a pp I cos 9, where 9 is 
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the relative angle between the normal vector to the soap 
film and the shear direction. To relate the normal drag 
force to the actual drag force of the Continuum Model, 
we need to project v 1 - in the shear direction (see Fig. 



11 



ii)), resulting in |v| = Iv^lcos^ = |v app |(cos 9) 2 




Fig. 11 Two projections are necessary to relate the veloc- 
ity v app of a soap film segment in the VF model to the local 
velocity v of the Continuum Model: (i) projection of the av- 
erage velocity of a soap film segment v app in the direction of 
the normal to that segment, and (ii) projection of the normal 
velocity of the soap film segment v ± back in the direction of 
shear. 8 is the angle between the normal vector to the soap 
film segment and the shear direction. It is the magnitude of 
the vectors that is displayed in the figure. 

Finally, we consider how the orientation of the soap 
films in our foam might affect the drag force. We assume 
that the foam is isotropic and proceed to average over 
all possible values of 9. As < (cos 6) 2 >= 1/2, our final 
expression for the continuum drag force coefficient c<j is 



giving the (continuum) drag force the required dimen- 
sions of force per area. 

In Sec. |4j we showed how the viscous stress has a 
XV dependence using dimensional arguments (see eq. [7| . 
Using these arguments, but rather defining the strain 



rate as a locally changing quantity (see eq. 20), we see 
that 



c v = c vl 1 ~ a A a -^ 2 X a 



(25) 



where the Herschel-Bulkley exponent a and the dimen- 
sionlcss quantity c v are free parameters. 

Using all of the approximations calculated in this sec- 
tion, we proceed to solve eq. [TT] numerically, accepting 
solutions only if they obey the inequality given by eq. 
19 as done in Sec. [5| The key difference here is that 
localisation length is a function of the product XV. 

The upper bound for the Continuum Model predic- 
tion is formulated in terms of XV by inserting eq.|24|and 



eq. 25 into eq. 16 resulting in 



l+(XV) 



2ac v a 1 ~ a A a - 1 / 2 \ 1+ 



(xvy 



(26) 



(24) 



(l + a)c d 

The corresponding lower bound must be found numeri- 
cally. To simplify this calculation, we fix A and allow V 
to vary. 

A comparison of the VF and Continuum Model re- 
sults can be seen in Fig. [12] where values of a = 0.3 
and c v = 0.26 are chosen as they give a reasonable 
prediction for both upper and lower bounds (although 
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0.2 < a < 0.4 gives a reasonable fit to the upper bound). 
The shaded region between these bounds indicates the 
range of all allowable localisation lengths, as predicted 
by the Continuum Model. Filled and open symbols rep- 
resent the VF simulation results, which are the same as 
in Fig. [7J only with the minimum localisation length of 
Irnin = A 1 / 2 subtracted to coincide with the Continuum 
Model predictions which give I = for V — ¥ oo. 



— upper, lower bounds 
• Sample 1 
A Sample 2 
■ Sample 3 

-4 Sample 4 
Sample 5 



natively, one may insert eq. [24] and eq. [25] into eq. [22] 
This gives 




0.05\ 0.1 0.15 0.2 

(XV) c = 0.056 M [y/J ^j 



Fig. 12 A comparison of the VF simulation results (open 
and filled symbols) and the prediction for the range of al- 
lowed localisation lengths as given by the Continuum Model 
(shaded region). A Herschel-Bulkley exponent of a = 0.3 and 
a stress overshoot of 20 % is found to give a good fit to the 
data. The critical cross-over point, (\V) C , as given by eq. 
|27| indicates the point below which the system yields a wide 
range of localisation lengths. 

The definition for the critical cross-over point, as 
given by eq. [22] may also be formulated in terms of Ay. 
This is achieved by inserting eq. [25] into eq. [21] and find- 



ing the point at which this line intersects eq. 26 Alter- 



V (l + a)c rf 

which is illustrated by the dotted lines in Fig. [12] The 
calculated critical cross-over point (AV") C = 0.056 7/ii 1 / 2 
fulfills its promise of offering a reasonable estimate of 
the point below which the system yields a wide range of 
localisation lengths. 

While the comparison between the VF Model and the 
Continuum Model presented in this section gives a fasci- 
nating theoretical explanation for the simulation results 
discussed, its details are far from precise. The location 
of the upper bound in Fig. [12] is simply an estimate, 
and further simulations may be needed to determine its 
exact location. In addition, many approximations were 
taken in relating the parameters of the two models. How- 
ever, it is remarkable that despite these approximations, 
a robust prediction can still me made. 

7 Outlook 

The apparent agreement of the simulation results in this 
paper with published experimental work suggests that 
the 2D VF Model may have further potential for describ- 
ing realistic foam dynamics. For more detailed studies to 
be conducted, however, the VF algorithm will need to be 
improved to decrease the required computation time for 
these types of simulations. Issues that we will address in- 
clude the effect of [ii {A) on localisation with the 2D VF 
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Model and on the value of the HB exponent. In addition, 
the dependence of the magnitude of the stress overshoot 
A on H2(A) will be investigated as it is critical to our 
understanding of its role as a mechanism for shear local- 
isation. It will be of interest to observe what happens to 
the location of the shear-band for samples with higher 
/12(A), in light of the results published in [33]. We also 
intend to compare our VF simulations with simulations 
using the Soft Disk Model [24]. 
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